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Abstract 

We demonstrate applications of algebraic techniques that optimize and certify polynomial inequal¬ 
ities to problems of interest in the operations research and transportation engineering communities. 

Three problems are considered: (i) wireless coverage of targeted geographical regions with guaran¬ 
teed signal quality and minimum transmission power, (ii) computing real-time certificates of collision 
avoidance for a simple model of an unmanned vehicle (UV) navigating through a cluttered environ¬ 
ment, and (iii) designing a nonlinear hovering controller for a quadrotor UV, which has recently been 
used for load transportation. On our smaller-scale applications, we apply the sum of squares (SOS) 
relaxation and solve the underlying problems with semidefinite programming. On the larger-scale 
or real-time applications, we use our recently introduced “SDSOS Optimization” techniques which 
result in second order cone programs. To the best of our knowledge, this is the first study of real-time 
applications of sum of squares techniques in optimization and control. No knowledge in dynamics 
and control is assumed from the reader. 

1 Introduction 

In this paper we consider applications of polynomial optimization in the area of operations research and 
transportation engineering. While techniques in more established areas of optimization theory such as 
linear, integer, combinatorial, and dynamic programming have found wide applications in these areas 
dH [ig [91 |22], the relatively newer field of polynomial optimization, which has gone through rapid 
advancements in recent years, may yet prove to reveal many unexplored applications. It is our aim in 
this paper to bring a few such applications to the attention of the operation research community and to 
highlight some algorithmic tools based on algebraic techniques that we believe are particularly suited 
for approaching problems of this sort. 

The fundamental problem underlying all of our applications is that of optimizing over nonnegative 
polynomials. This is the task of finding the coefficients Ca := Cai^...^an of some multivariate polynomial 
p{x) := p(xi,..., Xji) — in order to get p{x) > 0, either globally (i.e., Mx G M’^), or on certain 

basic semialgebraic sets. A basie semialgebraie set is a subset of the Euclidean space defined by a finite 
number of polynomial (in)equalities. That is a set of the form 

5 := {x G gi{x) > 0, hi{x) = 0}, 

where the functions hi are all multivariate polynomials. The polynomial optimization problem (POP) 
is itself a problem of this form. Indeed, the task of finding the minimum of a polynomial function q 
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on a basic semialgebraic set S is the same as that of finding the largest constant 7 such that q{x) — 7 
is nonnegative on S. There are, however, many other applications of optimization over nonnegative 
polynomials, some to be seen in this work. 

Our paper is organized as follows. In Section we briefly review the concept of sum of squares 
(sos) decomposition and its relation to semidefinite programming (SDP). This is a popular approach for 
certifying polynomial nonnegativity. While remarkably powerful, it often faces scalability limitations on 
larger-scale problems. As a potential remedy, we have recently introduced mm the concepts of diago¬ 
nally dominant and scaled diagonally dominant sum of squares (dsos and sdsos) decomposition, which 
instead of SDP result in linear programs (LP) and second order cone programs (SOCP) respectively. 
These concepts are also presented in Section 

In Section we consider the problem of providing guaranteed wireless coverage to certain basic 
semialgebraic subsets of the Euclidean space with minimum transmission power. The general problem 
here has been previously considered in the literature but we show that tools from polynomial optimiza¬ 
tion allow us to handle the problem in broader and arguably more realistic scenarios. Our next two 
examples are related to transportation problems in operations research. In particular, in Section we 
consider a simple model of an unmanned aerial vehicle (UAV), which aims to fly through a cluttered 
environment in a collision free manner. The techniques presented in this section can also be adapted for 
applications to ground vehicles. We demonstrate how one can choose a control law and at the same time 
find a formal certificate —an independently verifiable proof—that the resulting dynamics will guarantee 
no collisions with obstacles. We show that using our SOCP techniques, the underlying computational 
task can be carried out in the order of 20-30 milliseconds, hence making a plausibility claim about a 
real-time application of this approach. In Section we use the same technical tools to design a stabi¬ 
lizing controller for a quadrotor system, a device that has increasing potential for use in transportation 
(see Section]^. The designed controller prevents the quadrotor from losing balance when it is subject to 
environmental disturbance, or an external perturbation. The SDP resulting from this controller design 
problem is so large that it cannot be solved on our machine (3.4 GHz PC with 4 cores and 16 GB 
RAM). This is another example demonstrating the promise of our new sdsos machinery for handling 
problems of large scale. 

Our second and third applications include the employment of Lyapunov techniques to convert a 
problem in dynamics and control to a problem in polynomial optimization. Since we do not want to 
assume this background from the reader, we present the essentials of these very basic concepts in Section 
The mathematical background in this section (just like Section]^ is presented at a minimal level to 
make the paper self-contained, while keeping the focus on the applications and the algorithmic aspects. 
We end the paper with some brief concluding remarks in Section 

2 Algebraic certificates of nonnegativity via convex optimization 

The task of optimizing over nonnegative polynomials or even checking nonnegativity of a given polyno¬ 
mial, either globally or on a basic semialgebraic set, is known to be NP-hard [33]. This is true already 
for checking global nonnegativity of a quartic (degree-4) polynomial, or for checking nonnegativity of 
a quadratic polynomial on a set defined by linear inequalities. A popular relaxation scheme for this 
problem is through the machinery of the so-called sum of squares optimization. 

We say that a polynomial p is a sum of squares (sos), if it can be written as p = 'Yhi for some 
other polynomials qi. Obviously, such a decomposition is a sufficient (but in general not necessary [T9] ) 
condition for (global) nonnegativity of p. The situation where p is only constrained to be nonnegative 
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on a certain basic semialgebraic selQ 

5 := {x G M""! gi{x) > 0,.. .,gm{x) > 0} 

can also be handled with the help of appropriate sum of squares multipliers. For example, if we succeed 
in finding sos polynomials sq, si,..., such that 


p{x) = so(a:) + ^ Si{x)gi{x), (1) 

1=1 

then we have found a certificate of nonnegativity of p on the set S. Indeed, if we evaluate the above 
expression at any x G 5, nonnegativity of the polynomials sq, si..., 5 ^ imply that p{x) > 0. A Posi- 
tivstellensatz theorem from real algebraic geometry due to Putinar m states that if the set S satisfies 
the so-called Archimedean property, a property only slightly stronger than compactnes^ then every 
polynomial positive on S has a representation of the type Q, for some sos polynomials sq, 5 i,..., 5 ^ of 
high enough degree (see also [3l] for degree bounds). Even with absolutely no qualifications about the 
set 5, there are other Positivstellensatz theorems (e.g., due to Stengle EH) that certify nonnegativity 
of a polynomial on a basic semialgebraic set using sos polynomials. These certificates are only slightly 
more complicated than 0 and involve sos multipliers associated with products among polynomials gi 
that define S [36]. A great reference for the interested reader is the survey paper by Laurent [21]. 

The computational advantage of a certificate of (global or local) nonnegativity via sum of squares 
polynomials is that it can be automatically found by semidefinite programming. What establishes the 
link between sos polynomials and SDP is the following well-known theorem. Recall that a symmetric 
nxn matrix A is positive semidefinite (psd) if Ax > 0, Vx G M^, and that semidefinite programming 
is the problem of optimizing over psd matrices subject to affine inequalities on their entries [l3]. We 
denote the positive semidefiniteness of a matrix A with the standard notation A ^ 0. 

Theorem 2.1 (see, e.g., [35],[36]). A multivariate polynomial p{x) in n variables and of degree 2d is 
a sum of squares if and only if there exists a symmetric matrix Q (often called the Gram matrix) such 
that 

p{x) = z^Qz, . 

Q h 0, ^ ’ 

where z is the vector of monomials of degree up to d 

Z = [l,Xi,X2, . . .,Xn,XiX2, ■ ■ 


The search for the matrix Q satisfying a positive semidefiniteness constraint, as well as linear equality 
constraints coming from Q is a semidefinite programming problem. The size of the matrix Q in this 
theorem is 


n-\- d\ 


fn-\- d\ 

V d 


which approximately equals x While this number is polynomial in n for fixed d, it can grow 
rather quickly even for low degree polynomials. For example, the polynomials that we will be requiring 
to be SOS in our controller design problem for the quadrotor (Section]^ have 16 variables and degree 6 
and result in a Gram matrices with about half a million decision variables. A semidefinite constraint of 
this size is quite expensive—for example, the SDP solvers of SeDuMi [12] and MOSEK [2] fail to solve 
the quadrotor problem on our machine and quickly run out of memory. 


^In this formulation, we have avoided equality constraints for simplicity. Obviously, there is no loss of generality in 
doing this as an equality constraint h{x) = 0 can be imposed by the pair of inequality constraints h{x) > 0, —h{x) > 0. 

^In particular, if we have as an outer estimate a ball of some radius R in which our set S lives, then we can add a single 
quadratic inequality Xi < R to the description of S to have it satisfy the Archimedean property without changing the 
set. 
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2.1 DSOS and SDSOS Optimization 

In order to address the problem of scalability posed by SDP, we have recently introduced 1311 ! alterna¬ 
tives to SOS programming that lead to linear programs (LPs) and second order cone programs (SOCPs). 
The key insight there is to replace the condition that the Gram matrix Q be positive semidefinite (psd) 
with stronger sufficient conditions in order to obtain inner approximations to the cone SOSn^d of sos 
polynomials in n variables and of degree d. In particular, Q will be required to be either diagonally 
dominant (dd) or scaled diagonally dominant (sdd). We recall these definitions below. 


Definition 2.2. A symmetric matrix A is diagonally dominant fddj if an > '^j^i \ciij\ for all i. 
We will refer to the set of n x n dd matrices as DDn. 


Remark 2.1. It is clear from Definition |2.2| that the set DD^ has a polytopic description and can thus 
be optimized over using LP. 


Definition 2.3. Denote the set ofnxn symmetric matrices as S'^. Let ^ 2^2 ^ denote the symmetric 
matrix with all entries zero except the elements Mn^Mij^Mji^Mjj. Then, a symmetric matrix A is 
scaled diagonally dominant fsddj if it can be expressed in the following form: 


A — ^ 


Ain 

Mji 


Mij 

^33. 


h 0 . 


Remark 2 . 2 . The relationship between dd and sdd matrices is made clear in [5]. As we show there, a 
symmetric matrix A is sdd if and only if there exists a positive diagonal matrix D such that AD (or 
equivalently, DAD) is diagonally dominant. 

The set of n x n sdd matrices will be denoted by SDDn. We note that sdd matrices are sometimes 
referred to as generalized diagonally dominant matrices m- 

Theorem 2.4. The set of matrices SDDn can be optimized over using second order cone programming. 

Proof. Positive semidefiniteness of the 2x2 matrices in Definition |2.3| is equivalent to the diagonal 
elements Mn^Mjj, along with the determinant MnMjj — Mf^ being nonnegative. This is a rotated 
quadratic cone constraint and can be imposed using SOCP [ 6 ]. □ 

Remark 2.3. The fact that diagonal dominance is a sufficient condition for positive semidefiniteness 
follows directly from Gershgorin’s circle theorem. The fact that sdd implies psd is immediate from 
Definition |2.3| since a sdd matrix is a sum of psd matrices. Hence, denoting the set of n x n symmetric 
positive semidefinite matrices (psd) as we have from the definitions above that: 

DDn C SDDn C S+. 


We now introduce some naturally motivated cones that are inner approximations of the cone of 
nonnegative polynomials and that lend themselves to LP and SOCP. In analogy with the representation 
of SOS polynomials in terms of psd matrices (Theorem 2.1), we define the dsos and sdsos polynomials 
in terms of dd and sdd matrices respectively. 

Definition 2.5 ([5l[l]). 


• A polynomial p of degree 2d is diagonally-dominant-sum-of-squares (dsos) if it admits a represen¬ 
tation as p{x) = z'^{x)Qz{x), where z{x) is the standard monomial vector of degree d, and Q is a 
dd matrix. 
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• A polynomial p of degree 2d is scaled-diagonally-dominant-sum-of-squares (sdsos) if it admits a 
representation as p{x) — {x)Qz{x), where z{x) is the standard monomial veetor of degree d, 

and Q is a sdd matrix. 

We denote the set of polynomials in n variables and degree d that are dsos and sdsos by DSOS^^ 
and SDSOSn^d respeetively. 

The following inclusion relations are straightforward: 

DSOSn4 ^ SDSOSn4 ^ SOSn^^ 


Theorem 2.6. The set DSOSn^ is polyhedral and the set SDSOSn^ has a seeond order eone rep¬ 
resentation. For any fixed d, optimization over DSOS^^ (resp. SDSOSn^) can be done with linear 
programming (resp. seeond order eone programming), of size polynomial in n. 



Remark 2.4. While here we have chosen to define the DSOSn^ and SDSOSn^ cones directly in terms 
of dd and sdd matrices in order to expose their LP and SOCP characterizations, it is more natural to 
define them as sos polynomials of a partieular form. This alternate characterization is provided in [5]. 
In particular, we have the following equivalent definitions: 

• A polynomial p is dsos if it can be written as 

P = ^ Oiifnj + Pi, {mi + mjf + {mi - mjf, 

i hj 

for some monomials mi,mj and some constants > 0 . 

• A polynomial p is sdsos if it can be written as 

+ Y- 7/ rrijf, 

i hj 

for some monomials mi,mj and some constants 7 +, 7 ^“ > 0 . 

We will refer to optimization problems with a linear objective posed over the cones DSOSn^d^ 
SDSOSn^^ and SOSn^ as DSOS programs, SDSOS programs, and SOS programs respectively. In 
general, quality of approximation decreases, while scalability increases, as we go from SOS to SDSOS to 
DSOS programs. Depending on the size of the application at hand, one may choose one approach over the 
other. In this paper, we will be using SOS optimization (Section]^ and SDSOS optimization (Sections]^ 
and 1^ in our numerical experiments. The reader is referred to [2giai5] for many numerical examples 
involving DSOS optimization. We also remark in passing that SDSOS or even DSOS programming enjoy 
many of the same theoretical (asymptotic) guarantees of SOS programming—results of this nature are 
proven in [5]. 

We now proceed to some potential operations research applications of the tools discussed so far. 
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3 Wireless coverage with minimum transmission 

In the problem considered in this section we have a number n of wireless electromagnetic transmitters 
located at positions on the plane. Each transmitter is an omnidirectional power 

source, emitting waves in all directions with equal intensity. Due to the laws of electromagnetics, the 
energy Ei propagated from each jamming device is inversely proportional to the squared distance from 
the device: 

Et{x,y) j-y _-.^2 ’ 

where A is some propagation constant, set hereafter to 1 with no loss of generality, and q is the 
transmission rate of device i. The goal is to make sure that certain regions of the plane are guaranteed to 
receive a given cumulative energy level of at least C units, while minimizing transmission power. These 
regions can for example be populated urban geographical domains where a wireless service provider 
would like to guarantee a certain level of signal quality. 

The problem we describe is motivated by some interesting and relatively recent work in US], US] (see 
also the thesis jl4|). where the motivation is instead to jam the communication network of an adversary 
with a wireless transmitter. We note, however, that there are a few differences between our setting and 
that of m and [15], the main one being the assumption about the region to be covered. Reference m 
assumes that this region is a set of isolated points (the location of the adversary is known) and this 
results in a simplified problem. However, more complex objectives are considered by the authors; e.g., 
the goal is to make the communication graph of the enemy disconnected, or to jam a prescribed fraction 
of the enemy locations, or to decide which transmitters to turn off. On the opposite end, the work 
in m assumes absolutely nothing about the location of the adversary. As a result, the goal is to cover 
an entire rectangular region by a prescribed level of jamming power. Our setting, by contrast, allows 
for the region to be covered to be the union of arbitrary basic semialgebraic sets (see, e.g.. Figure [T(a)] ); 
this obviously enhances the modeling power. We should also comment that neither our work, nor the 
works in m and [15], satisfactorily address the more difficult problem of optimizing over the location 
of the transmitters. 

A formal summary of our setting is as follows. We are given as input the following quantities: C 
(required coverage level), 7 ^ (upper bounds on transmission rates), (x^, ^^), i = 1 ,..., n (location of our 
transmitters), Bj^j = 1,... ,m (basic semialgebraic sets describing regions to be covered). We assume 
that the transmitters are outside of the location sets Bj. The goal is to find transmission rates q to 
solve the following optimization problenj^ 

minimize Y17=i 

Ci < -fi, Vi = l,...,n, (3) 

E{x,y) := YJl=i > C, V(x, y) G j = 1 ,..., m. 

Note that the latter constraints are requiring certain rational functions to be nonnegative on certain 
basic semialgebraic sets. Upon taking common denominators, we can rewrite these constraints as 
polynomial inequality constraints: 

p(x,y) := - Xif + (9 - 5,)2| + Y.U c. X[kd(x - %)" + is - Skf] > 0, , 

V(x,y) € Bj,j ^ 1,... ,m. 

Observe that the degree of the polynomial p(x, y) is two times the number of transmitters. Since we 
are dealing with polynomial inequalities in only two variables, we have no scalability issues restraining 
us from applying the sos relaxation. Let each set Bj be defined as 

Ej = {a;| ^ > 0}, 

®An alternative reasonable objective is to minimize ma:xi a. This can as easily be handled. 


6 








(a) The input to the problem. 



(b) Log of the cumulative energy at (c) Region that receives a total en- 
the SDP solution. ergy of at least C = 10 units (in 

dark red). 


Figure 1: An instance of the wireless coverage problem, 
for some bivariate polynomials The optimization problem that we will be solving is: 


minimize 


ELi 

k • 

= ^0 + '^iLl ^j,k9j,i 

where p is as in Q and (JQ^aj^k bivariate polynomials whose degree is upper bounded by some 
even integer d. Note that the above is a semidefinite programming problem (via Theorem 2.1) with 
decision variables consisting of the scalars q and the coefficients of the polynomials ao^aj^k- If is easy 
to see that for each value of the degree d, the optimal value of (§ is an upper bound on the optimal 
value of Moreover, since in our setting each set Bi satisfies the Archimedean propert}Q Putinar’s 
Positivstellensatz tells us that by increasing d, we will be able to solve (§ to global optimality. 

Let us now solve a concrete example. Our input data is demonstrated in Figure [T(a^ We have two 
transmitters, located at points (1,1.5) (called transmitter 1) and (2,1) (called transmitter 2) on the 
plane. The area to be covered is given by the five ellipsoidal regions 


P 


j = l,...,m, 

SOS, 


(5) 


Bj = {z := {x,yf\ {z - ZjfAj{z - Zj) < 


with Ai = 


3 1' 
_1 1_ 

5^2 = 

'o CO ' 

T—1 O 

1_1 

5^3 = 

' O 

o 

1_1 

II 

■ 1 -1' 
-1 3 _ 

II 

' O CJi' 

o 

1_1 


,zi = (l.l,1.75y,Z2 = 


(1.25, 2)’^, zs = (1.5,1.75)"^, 2:4 = (1.8,1.8)"^, 25 = (2,1.4)’^, ai — a 2 = = a 4 — 0.1,05 = 0.2. 

The required energy level on these areas is C = 10 and the upper bounds on both transmission 
rates ci, C 2 is 11. We first would like to know if by only turning on one of the two transmitters we can 
meet the required energy level. For the transmitter at the location (2,1), the optimal value of the SDP 
in ^ with degree of sos multipliers set to zero (i.e., constant multipliers) is 17.594. In fact, in this 
case, we know that this upper bound is already exact! This is because in the case of one transmitter, 
the polynomial p in (|^ is quadratic. If a quadratic polynomial is nonnegative on a region defined by 
another quadratic, this fact is always certified by a constant degree multiplier—this is the celebrated 
5-lemma; see m- Similarly, if we solve the problem for the transmitter located at (1,1.5), the optimal 
value of (§ which matches the optimal value of (© is 11.446. So our task is indeed not achievable with 
one transmitter only. 

With both transmitters on, the SDP in ^ is infeasible for degree-0 sos multipliers (giving an upper 
bound of infinity). However, when we increase the degree of these multipliers to 2, a solution is returned 
with Cl = 2.561 and C 2 = 5.550 at optimality. By further increasing the degree of our sos multipliers. 


^Indeed each set Bi is compact and the entire environment can be placed in a ball of some prescribed radius R. This 
quadratic constraint can be added to the description of each Bi to satisfy the Archimedean property. 
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no improvement in optimal value is observed and we conjecture that the numbers above are already 
optimal for the original problem Q. Figure [l(b)| shows the logarithm of the cumulative energy level 
E{x^y) at each point in space. (The logarithm is taken to better observe the dispersion of energy.) 
Figure [T(^ shows all pixels that receive the required energy level of C = 10 units. As promised, all five 
ellipsoids are covered and interestingly the boundary of the region covered touches two of the ellipsoids. 

4 Lyapunov theory and optimization 


,- Lyapunov 

' Questions about 1 Theory ^ SOS 

' . . I X « _ Search for sign- ' ■ Polynomial I ■ Semidefinite ' 

I , ^ ^ , r definite functions . inequa ities : programming ' 

I (e.g., stability, safety) 1 l -> i_ 2-^ __ tj 

^SDSOS 

! socp 1 

I — — — — _ J 

Figure 2: The steps involved in Lyapunov analysis of dynamical systems via convex optimization. 

The examples presented in our next two sections involve decision-making about trajectories of dy¬ 
namical systems. The machinery that allows us to reduce such tasks to problems in optimization is 
Lyapunov theory. As depicted in Figure the general idea is the following: In order to guarantee that 
trajectories of dynamical systems satisfy certain desired properties, it will be enough to find certain 
scalar valued functions that satisfy certain inequalities. These functions will be parameterized as poly¬ 
nomials and DSOS/SDSOS/SOS relaxation techniques will be used to find their unknown coefficients 
in such a way that the desired inequalities are automatically satisfied. For our two applications, we 
explain next what these inequalities actually are. 

Barrier functions (Section [^. Consider a differential equation x = /(x), where x denotes the 
derivative of the state vector x with respect to time and / : ^ is a polynomial function. Suppose 

we are given two basic semialgebraic sets and 5^nsafe want to guarantee that trajectories 

starting in would never end up in 5^nsafe‘ This guarantee can be achieved if we succeed in finding 

a function V : ^ M, called a harrier funetion ISH!, m, 0, with the following three properties: 

V{x) <1 \/xe 5gafe, V >1 Vx G < 0 Vx. 

The expression V denotes the time derivative of V along trajectories. If V is a polynomial, V will also 
be a polynomial given (via the chain rule) by: 

V{x) = {VV{x)J{x)). 

The three inequalities above imply that it is impossible for a trajectory to go from to 5^nsafe 

since the function V evaluated on this trajectory would need to go from a value less than one to a value 
more than one, but that cannot happen since the value of V is non-increasing along trajectories. 

Stability and region of attraction computation (Section]^. Suppose once again that we have 
a differential equation x = f{x) with origin as an equilibrium point (i.e., satisfying /(O) = 0). In 
numerous applications in control and robotics, one would like to make sure that deviations from an 
equilibrium point tend back to the equilibrium point. This is the notion of asymptotic stability. A 
particularly important problem in this area is the so-called “region of attraction (ROA) problem”: For 




















what set of initial conditions in do trajectories flow to the origin? This question can be addressed 
with Lyapunov theory. In fact, Lyapunov’s stability theorem (see, e.g., [23l Chap. 4]) tells us that if 
we can find a (Lyapunov) function V : ^ M, which together with its gradient VV satisfies 


y(x)>0 Vx ^ 0, and V{x) = {W{x)^ f{x)) < 0 Vx G {x| y(x) </?, x ^ 0}, (6) 

then the sublevel set {x\ V{x) < 13} is part of the region of attraction. Notice again that if / is a 
polynomial function (an immensely important case in applications [3l Chap. 4]), and if we parameterize 
y as a polynomial function, then the search for the coefficients of V satisfying the conditions in Q is 
an optimization problem over the set of nonnegative polynomials. 

5 Real-time Planning with Barrier Functions 



Figure 3: An illustration of the states of the UAV model we consider. 

One promising application domain for polynomial optimization in transportation is for real-time plan¬ 
ning and control on autonomous vehicles. In this example, we consider such an application for a simple 
model of an unmanned aerial vehicle (UAV) navigating through a cluttered two dimensional environ¬ 
ment. In order to make the navigation task more realistic, we also consider a bounded but uncertain 
“cross-wind” term in the dynamics. This results in an uncertain differential equation and requires rea¬ 
soning about families of trajectories that the system could end up following, making the problem more 
challenging. The states and dynamics of the UAV are inspired by the widely-used Dubins car model 
m and are given by: 


X 


X 


—V sin?/ -h w 

y 

, x = f(x,u,w) = 

y 

= 

V cos 

. ^ . 


. ^ . 


u 


where x and y are the x and y positions of the UAV in the environment, v = 1 m/s is the speed of 
the airplane, is the yaw angle, u is the control input and w is the “cross-wind” (bounded between 
[—0.05, 0.05]). An illustration of the states of the model are given in Figure]^ We Taylor expand these 
dynamics to degree 3 to obtain polynomial dynamics in order to use DSOS/SDSOS/SOS programming. 



Figure 4: A barrier function computed for a particular initial state and obstacle configuration. The 
UAV is guaranteed to remain safe when the controller is executed despite the effects of the cross-wind. 
The green curve is a level set of a degree-4 polynomial found by SDSOS optimization. 
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Figure 5: UAV successfully navigating through different obstacle environments using the planning algo¬ 
rithm described in Section]^ A video of the navigation can be found at http://youtu.be/J3a6v0tlsD4, 
This video also shows the barrier certificates (not shown here) as they get updated in real time. 

Our goal is to make the UAV navigate through cluttered environments that are unknown pre-runtime 
while avoiding collisions with the obstacles in the environment despite the effects of the cross-wind on 
the vehicle dynamics. In order to achieve this, we pre-compute five control primitives that the UAV 
can choose from at runtime. These controllers take the form: 

Mi(x) = -K{'lp- i>des,i), i = 1,..., 5. (8) 

These control primitives cause the UAV to control its yaw angle to a particular angle We 

choose iF = 50 and ^l^des.i — 0 rad, il^des ,2 = —207r/180 rad, ipdes.s = 207r/180 rad, ipdes,4. = —457r/180 
rad, ipdes,5 = 457r/180 rad. 

The UAV’s task is to choose from these control primitives in order to navigate its way through the 
environment. After executing a particular chosen primitive for a short interval of time (1/20 seconds 
in our case), the UAV replans by choosing a control primitive again. Hence, the key decision that our 
planning algorithm needs to make is the choice of control primitive given a particular configuration 
of obstacles in its environment. We take our inspiration from [8], which uses barrier eertifieates for 
verifying the safety of a controller given a particular set of obstacles (but does not consider the case 
where obstacle positions are not known beforehand and decisions must be made in real time). Similarly, 
other previous SOS programming approaches m to collision avoidance have involved solving SOS 
programs offline and then using these precomputed results to do planning in real time. In contrast, in 
our example here, the optimization problems are solved in real-time. We describe our approach below. 

At every control iteration, we identify the closest two obstacles in the environment in front of the 
UAV. We then evaluate each control primitive Ui and check if executing it from our current state will 
result in the UAV avoiding collision with the obstacles. The first safe controller found is executed. The 
safety of a controller can be checked by computing barrier functions using the polynomial optimization 
approaches described in Section]^ Denoting the current state as xq = (xo^yo^^o) and the obstacle sets 
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as Xoi)s^i C and Xobs,2 C we use polynomial optimization to search for a function V{'x.) of degree 
4 that satisfies the following conditions: 


y(xo) = 0, 


(9) 

V^(x) > 1, 

^ Xobs,i’) ^ 

(10) 

U(x,ie) = 

dV 

<0, Vx G X, Vie G [—0.05,0.05]. 

C/X 

(11) 


Here, X is a “large” set that the system is guaranteed to remain within for the duration of time for which 
the control primitive is executed. In particular, we choose it to be the unit sphere around the current 
state. The conditions above imply that the state x is constrained to evolve within the 1-sublevel set 
(in fact the 0-sublevel set) of the function V (x) and is thus guaranteed to not collide with the obstacles 
despite the effects of the cross-wind. 

Hence, at each control iteration we need to solve a maximum of 5 optimization problems, all of which 
are independent and can be parallelized. In our example, we use SDSOS programming to compute 
barrier functions and observe running times of approximately 0.02 — 0.03 seconds for feasible problems 
and 0.08 — 0.09 seconds for infeasible problems (i.e., problems where no barrier function can be found) 
using the Gurobi SOCP solver [1] (a more thorough running time analysis is presented later). Hence, 
a real-time implementation of this approach on a hardware platform is plausible. Such a hardware 
implementation can benefit from already-existing SOCP solvers that are specifically designed to run 
on embedded systems m, EH]. In particular, m presents an approach for generating stand-alone C 
code for an SOCP solver that can run very efficiently and with low memory footprint. The use of such 
real-time SOCP solvers has already been considered for tasks such as landing of spacecraft (e.g., for 
NASA’s Mars exploration project) [12]. 

A particular example of a barrier function computed for the controller ui is shown in Figure 
The obstacles are shown in red and the initial state of the UAV is also plotted. The 1-level set of the 
computed barrier is plotted in green and certifies that the initial state is guaranteed to remain safe when 
the controller is executed. 

Figure [^demonstrates the performance of the algorithm described above with SDSOS programming 
used to compute barrier functions on a number of environments. Each subfigure shows a randomly cho¬ 
sen environment (with obstacle positions chosen from the uniform distribution) with circular obstacles 
that the UAV has to navigate. The trajectory traversed by the UAV following the described planning 
algorithm is indicated in these plots and remains collision free in each case. Note that the original 
(non-Taylor expanded) dynamics are used for the simulations. 

We end the discussion of this example by comparing running times and performance of the SDSOS 
and SOS approaches to this problem. In order to do this, we fix the initial state of the vehicle to 
be (O,O,'0o) foi* varying values of t/^o- For each we randomly sample 100 different environments 
containing two obstacles each. The obstacles are disks of radius 0.03 m with centers (xc, yc) uniformly 
sampled in the range Xc G [—0.2, 0.2] m, yc G [0,0.2] m. For each environment, we attempt to find 
a valid barrier certificate for the first controller in our library (i.e., the one that servos the vehicle to 
"^des,! = 0)- The results are summarized in Table [^ which presents the number of environments (out of 
100) for which a barrier certificate was successfully found using SDSOS and SOS programming. As the 
table illustrates, the number of times SDSOS programming fails to find a barrier certificate when SOS 
programming sueeeeds is quite small. 

We also compare running times of the two approaches in Figure We use the Gurobi SOCP solver 
[1] for the SDSOS problems and SeDuMi [12] as the SDP solver for SOS problems. As the histograms 
of running times illustrate, the SDSOS approach is significantly faster than the SOS approach. We note 
that while the MOSEK SOCP/SDP solvers are typically faster, we were unable to make these work on 
this problem due to numerical issues. 
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Table 1: Comparison of percentage of times a valid barrier certificate was found using SDSOS and 
SOS programming for randomly sampled obstacle environments and initial yaw angles. (Only the ratio 
between the two is meaningful here.) 



running time (s) 


Figure 6: Histograms of running times for SOS and SDSOS approaches on the collision avoidance 
problem. 


6 Nonlinear Control Design for a Quadrotor Model 

Quadrotors (see Figure have recently been recognized as a popular platform for academic research 
in systems theory due to their agile maneuvering capabilities and inexpensive cost 1311 EOj. They have 
also been considered for the task of load transportation, not only in laboratory settings (320 but also 
by the aerospace companies Bell and Boeing and the online retail company Amazon^ In this section, 
we consider the problem of designing a nonlinear stabilizing feedback controller for the quadrotor’s 
hovering configuration, which is relevant to almost all of its applications. In addition to a stabilizing 
controller, we also obtain a formal certificate of stability of the resulting system. This certificate takes 
the form of an inner approximation of the region of attraction (ROA), i.e., the set of initial conditions 
the controller is guaranteed to stabilize to the goal position. 



Figure 7: We design a hovering controller for the quadrotor model described in [30]. (Image from [30].) 

We use the dynamics model described in [30] for our numerical experiments. The model includes 16 

video corresponding to the paper is available at https://www.youtube.com/watch?v=YBsJwapanWI 
^ https://www.youtube.com/watch? v=Le46ERPMlWU 
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states: 


X := lxi,y,z,(fi,0,'ifi,x,y,z,p,g,r,u>i,u>2,u>3,u;4j, 

where xi^y^z are the coordinates of the center of mass of the system, are the Euler angles 

describing its orientation, g, r are angular velocities of the quadrotor expressed in the body frame, 
and uji^i = 1, ...,4, are the angular speed of the rotors. The rotor angular speeds cannot be controlled 
directly and have nontrivial dynamics. The control inputs of the system are thus the desired speed of 
the rotors (the rotors take some time to catch up to the desired speed). 

In the end our system takes the form x = f{x) -\- g{x)u{x) with / and g given and the control u as 
a decision function. We use the method presented in our earlier work [25] in collaboration with Russ 
Tedrake to design a hovering controller u for the system. The fixed point corresponding to the hovering 
configuration has the first twelve states of the system equaling 0 but with non-zero rotor speeds uji 
counteracting the force of gravity. The dynamics of the system are Taylor expanded to degree 3 in 
order to obtain polynomial dynamics. We search for a degree 2 Lyapunov function V{x) and a degree 3 
feedback controller u{x) in order to maximize the size of the region of attraction (ROA) of the resulting 
closed-loop system (i.e., the differential equation with u{x) plugged in). We use SDSOS programming 
since the state space is too large for SOS programming to handle, causing our computer to run out of 
memory. The resulting optimization problem is: 


max p (12) 

p,L{x),V {x),u{x) 

s.t. V{x) € SDSOSie ,2 

-V{x) + L{x){V{x) -p)e SDSOSi&^q 
L{ x) € SDSOSi&,4 
Y^V{eP = l. 

3 


Here, L{x) is a nonnegative multiplier term and ej is the j-th standard basis vector for the state 
space From our discussion in Section it is easy to see that the above conditions are sufficient 

for establishing Bp — {x ^ | V{x) < p} as an inner estimate of the region of attraction for the 

system. When x G Hp, the second constraint implies that R(x) < 0 (since L{x) is constrained to be 
nonnegative). The last constraint normalizes V{x) so that maximizing the level set value p leads to 
enlarging the volume of the ROA. 

The optimization problem (12) is not convex in general since it involves conditions that are bilinear 
in the decision variables. However, problems of this nature are common in the SOS programming 
literature (see e.g. m) and are typically solved by iteratively optimizing groups of decision variables. 
Each step in the iteration is then a SDSOS program. This iterative procedure is described in more detail 
in [25] and can be initialized with the Lyapunov function from a Linear Quadratic Regulator (LQR) 
controller [7]. The iterations are terminated when the objective changes by less than 1 percent. 



(a) X - y subspace 



(b) X - z subspace 



(c) X - pitch subspace 


Figure 8: Slices in different subspaces of the hovering ROA of the quadrotor system. 
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An important observation is that unlike the sets POSn^d SOSn^^ the sets DSOSn^d s^^nd 
SDSOSn 4 are not invariant to coordinate transformations, i.e., a polynomial p{Ax) is not necessarily 
dsos (resp. sdsos) even if p{x) is dsos (resp. sdsos). Thus, performing coordinate transformations on 
the problem data (e.g., on the state variables of a dynamical system) can sometimes have an important 
effect. We describe a particular coordinate transformation that is intuitive and straightforward to imple¬ 
ment. It can be used for problems involving the search for Lyapunov functions, and can potentially be 
extended to other problems as well. In particular, given a Lyapunov function V{x) we find an invertible 
affine transformation that simultaneously diagonalizes the Hessians of V{x) and —V{x) evaluated at the 
origin (this is always possible for two positive definite matrices). The intuition behind the coordinate 
change is that the functions V{x) and —V{x) locally resemble functions of the form x^Dx (with D 
diagonal), which are dsos polynomials that are “far away” from the boundary of the DSOS (and hence 
SDSOS) cone. We solve the optimization problem (12) after performing this coordinate transformation. 
The transformation is then inverted to obtain ROAs in the original coordinate frame. 



Figure 9: A sampling of five initial conditions that are stabilized by our controller. The goal position is 
shown in green, the stabilized initial conditions in red, and the intermediate trajectories in blue. 

Each iteration of the algorithm employed for solving the optimization problem ( [T^ takes approxi¬ 
mately 15 minutes, with convergence occurring between 15 and 20 iterations. Figure shows slices of 
the computed ROA in multiple subspaces of the state space. As the plot illustrates, we are able to verify 
stability of the closed loop system for a large set of initial conditions. A qualitative demonstration of 
the performance of the controller is given in Figure]^ The system is started off from five different initial 
conditions (shown in red) and our nonlinear hovering controller is applied. The resulting trajectory is 
shown in blue. In each case the quadrotor is able to stabilize itself to the goal configuration (green). 

7 Conclusions 

In this paper, we demonstrated three applications of optimization problems over the set of nonnegative 
polynomials that may be of interest in operations research and transportation engineering. We hope 
to have conveyed the message that the problem of certifying polynomial inequalities appears in more 
diverse areas than one might think. There are powerful tools for approaching this problem based on 
the sum of squares relaxation and semidefinite programming. We believe that our recently introduced 
techniques of DSOS and SDSOS optimization, which are LP and SOCP-based alternatives to sum of 
squares programming, can pave the way to new applications of algebraic techniques in optimization—in 
particular, applications that are large-scale or real-time. 
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